#include "aa_util.h"
#include "../pdb_index/aa_table.h"
#include <string.h>
#include <assert.h>
#include <math.h>
off_t find_in_aa_by_name(const struct store * store, const off_t start_pos, const char * name) {
  int i;

  struct atom * atom = store_read(store, start_pos);
  int start_residue = atom->residue;
  int last_residue = start_residue;
  // Needed because residues may be numbered zero, so we need to use
  // the atom serial number for the EOF check...
  int last_atom_serial_number = atom->serialNumber;  
  store_release(store, start_pos);

  
  /* Unfortunately, some pdb files (1a6r, for instance) have residues
     with negative indices.  This has the side effect of making
     things more complicated.... 

  // Very bad things could happen if the
  // start residue is zero...This function
  // might not terminate, and could grow the
  // store until the OS killed the program
  // off... Since unix has sparse files, you
  // could see a 'multi-terabyte' file within
  // seconds...better leave this assert here.

  assert (start_residue); */

  assert (last_atom_serial_number);

  
  for(i = start_pos; i >= 0 && last_residue == start_residue; i--) {
    atom = store_read(store, i);
    last_residue = atom->residue;

    // Don't set this. (don't use it for anything, and we want to
    // preserve the inital value for the next loop.

    //    last_atom_serial_number = atom->serial; 

    if(last_residue == start_residue && !strncmp(name, atom->name, 4)) {
      store_release(store, i);
      return i;
    }
   
    store_release(store, i);
  }

  


  last_residue = start_residue;
  //  last_atom_serial_number = start_serial_number // We don't need this because last_atom_serial_number already equals the start_serial_number...

  // Assuming that the file is tightly packed, when we run off the end
  // of it, last_atom_serial_number will default to zero.

  for(i = start_pos + 1; last_atom_serial_number && (last_residue == start_residue); i++) { 
    atom = store_read(store, i);
    last_residue = atom->residue;
    if(last_residue == start_residue && !strncmp(name, atom->name, 4)) {
      store_release(store, i);
      return i;
    }
    store_release(store, i);
  }
  return -1;

}

void * filter_atom_by_name(struct atom* a, char atomName[5]){

  if(strncmp(a->name, atomName, 4)) {
    return 0;
  } else {
    return a;
  }
}
/**
   _REMOVES_ atoms that are in the given residue..
*/
void * filter_atom_by_residue_number(struct atom* a, int residueNumber){
  if(a->residue == residueNumber) {
    return 0;
  } else {
    return a;
  }
}


void find_aa_range(const struct store * atomStore, const off_t i, off_t * start, off_t * stop) {

  struct atom * a = store_read(atomStore, i);
  int resNumber = a->residue;
  int thisRes = resNumber;
  store_release(atomStore, i);
  for(*start = i; (*start >= 1) && (thisRes == resNumber); (*start)--) {
    a = store_read(atomStore, *start-1);
    thisRes = a->residue;
    store_release(atomStore, *start-1);
  }
  *start = *start + 1;

  thisRes = resNumber;

  for(*stop = i; thisRes == resNumber; (*stop)++) {
    a = store_read(atomStore, *stop+1);
    thisRes = a->residue;
    store_release(atomStore, *stop+1);
  }
  *stop = *stop - 1;
  
}

int can_calculate_aa_mass_center(const struct store * atomStore, const off_t start, const off_t stop) {
  off_t i = start;
  int k1,k2,k3,k4,k5;
  while(i<stop) {
    struct atom * a = store_read(atomStore, i);
    k1 = abs(strncmp(a->name, " N  ", 4));
    //    k2 = abs(strncmp(a->name, " CA ", 4));
    k3 = abs(strncmp(a->name, " O  ", 4));
    k4 = abs(strncmp(a->name, " C  ", 4));
    k5 = abs(strncmp(a->element, " H",2));
    store_release(atomStore, i);
    if((k1)/*&&(k2)*/&&(k3)&&(k4)&&(k5)){
      return 1;
    }
    i++;
  }
  return 0;
}


void calculate_aa_mass_center(const struct store * atomStore,
			   const off_t start, const off_t stop,
			   struct xyz_int * output) {
  off_t i = start;
  double weight, sum_weight=0,x=0,y=0,z=0;
  int k1,/*k2,*/k3,k4,k5;
  
  while(i <= stop){
    struct atom * a = store_read(atomStore, i);
    k1 = abs(strncmp(a->name, " N  ", 4));
    //    k2 = abs(strncmp(a->name, " CA ", 4));
    k3 = abs(strncmp(a->name, " O  ", 4));
    k4 = abs(strncmp(a->name, " C  ", 4));
    k5 = abs(strncmp(a->element, " H",2));
    if((k1)&&/*(k2)&&*/(k3)&&(k4)&&(k5)){
       weight = find_atom_weight( a->element);
       x += (double)weight * (int)(a->x);
       y += (double)weight * (int)(a->y);
       z += (double)weight * (int)(a->z);
       sum_weight += weight;
    }
    store_release(atomStore, i);
    i++;
  }
  
  output->x =(int) (x/sum_weight +0.5);
  output->y =(int) (y/sum_weight +0.5);
  output->z =(int) (z/sum_weight +0.5);
}

#define OK 1
#define NONE 0
#define ERROR -1

int calculate_aa_functional_group(const struct store * atomStore, const off_t start, const off_t stop, struct xyz_int * out) {
  struct atom * a = store_read(atomStore, start);
  struct atom * b;
  int a_idx = aa_get_table_offset(a->residue_type);
  int ret;
  const char * name_1 = aa_active_sites[a_idx].first_atom;
  const char * name_2 = aa_active_sites[a_idx].second_atom;
  off_t atom1 = find_in_aa_by_name(atomStore, start, name_1);
  off_t atom2 = find_in_aa_by_name(atomStore, start, name_2);

  store_release(atomStore, start);

  a = 0;
  b = 0;

  if(atom1 != -1) {
    a = store_read(atomStore, atom1);
  }
  if(atom2 != -1) {
    b = store_read(atomStore, atom2);
  }
  
  if(name_2[0]) {
    if((!a) || (!b)) {
      ret = ERROR;
    } else {
      out->x = (int)((a->x + b->x) / 2.0);
      out->y = (int)((a->y + b->y) / 2.0);
      out->z = (int)((a->z + b->z) / 2.0);
      ret = OK;
    }
  } else if(name_1[0]){
    if(!a) {
      ret = ERROR;
    } else {
      out->x = a->x;
      out->y = a->y;
      out->z = a->z;
      ret = OK;
    }
  } else {
    ret = NONE;
    out->x = 0;
    out->y = 0;
    out->z = 0;

  }
  if(a) {
    store_release(atomStore, atom1);
  }
  if(b) {
    store_release(atomStore, atom2);
  }

  return ret;
  
}
